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SECTION  I 


INTRODUCTION 

One  of  the  ubiquitous  problems  in  the  assessment  of  fatigue  life  is 
the  determi lation  of  the  "strength"  reitaining  when  structural  components 
have  been  si  ejected  to  a  specified  stress  regime  for  a  known  length  of  time. 

Even  in  the  simplified  circumstances  when  nominally  identical  structural 
components  are  subjected  to  only  an  alternating  stress  of  a  given  maximum 
value  with  fixed  stress  variation  (such  as  might  be  the  case  for  a  helicop¬ 
ter  rotor)  the  fatigue  life  will  still  exhibit  sufficient  variability  that 
the  mathematical  mudel  deemed  most  suitable  is  that  the  fatigue  life  is  a 
random  variable  with  an  appropriate  distribution.  Of  course,  the  resolu¬ 
tion  of  the  question  of  which  parametric  class  of  distributions  is  the  most 
useful  and  appropriate  has  not  been  universally  agreed  upon. 

A  concept  which  has  historical  precedent  is  that  the  characteristic 
life  8  of  a  structural  component,  which  has  been  selected  by  chance  and 
then  subjected  to  repeated  applications  of  a  maximum  stress  level  s  , 
can  be  expressed  as  a  function  of  that  stress,  say 

6  -  K(s)  (1) 

where  K  is  a  decreasing  function,  usually  assumed  to  be  known  except  for 
some  parameters. 

The  equation  (Eq  (1))  Is  called  the  Wohler  equation  and  it  expresses  the 
characteristic  life  in  cycles  to  failure  as  a  function  of  the  maximum  stress. 
In  statistical  parlance  this  is  a  regression  equation.  The  plot  of  sample 
data  used  to  estimate  this  relation  is  called  the  statistical  S-N  diagram 
i.e.  stress  versus  number  of  cycles  to  failure.  Various  specifications  have 
been  adopted  for  K  .  One,  for  example,  is  the  Liindberg  formula 
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(2) 


K(s) 


( 
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s-s. 
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l. 
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for  s  >  s 

0 


where  is  the  value  of  the  stress  called  the  "fatigue  limit"  and  the 

constants  a,b  >  0  are  determined  from  the  material.  These  parameters  are 
usually  estimated  by  least  squares,  or  by  eye,  from  the  statistical  S-N 
diagram. 

One  approach  is  to  specify  a  life  length,  before  which  failure  will 
occur  with  preassigned  probability  P  according  to  the  postulated  distri¬ 
bution  of  life,  and  to  construct  a  P-S-N  diagram  (see  Bastenaire  [  1  ]  for 
notation  and  terminology)  in  accord  with  a  cumulative  damage  hypothesis, 
such  as  Miner's  Rule,  using  equi-probability  curves. 

Miner's  Rule  utilized  with  the  Wohler  regression  equation  would  be  as 
follows:  In  a  loading  spectrum  which  contains  n^  cycles  of  maximum  stress 

level  s^  for  i  *  l,...,k  the  characteristic  life  is  given  (in  number  of 
spectra  which  can  be  repeated) by  the  reciprocal  of 


k 


I 


n. 


K(si) 


This  brief  historical  perspective  is  presented  to  point  out  how  the 
current  approach  attempts  to  overcome,  or  strengthen,  the  weakness  of  the 
theory  previously  employed. 

The  first  point  is  that  the  approach  just  mentioned  does  not  recognize 
the  decided  influence  upon  fatigue  life  which  the  load  order  within  the 
spectrum  plays.  This  load  order  influence  has  been  known  to  research 
workers  fr r  a  long  time,  see  Ref  [2],  but  the  implementation  of  such  knowledge 
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into  practical  formulae  useable  in  preliminary  airframe  design  has  not  been 
accomplished. 

The  second  point  is  that  there  is  n  great  danger  in  assuming  a  theory 
which  has  as  a  necessary  component  a  so-called  "fatigue  limit"  i.e.  a  value 
below  which,  if  the  stress  is  reduced,  the  fatigue  life  becomes  infinite. 
Such  a  theory  does  not  take  into  account  flaws,  accidents  or  "the  thousand 
natural  shocks  that  metal  (and  flesh)  is  heir  to"  which  can  and  do  initiate 
fatigue  failure. 
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SECTION  II 


SOME  CONSIDERATIONS  IN  THE  DERIVATION  OF 
LIFE  LENGTH  MODELS  FOR  FATIGUE 

The  simple  model  which  was  originally  utilized  to  obtain  the  distribu¬ 
tions  derived  in  Ref  [3]  was  as  follows :  During  the  stable,  microscopic  phase 
of  fatigue  crack  growth,  following  crack  initiation,  the  incremental  crack 
extension  during  each  repetition  of  the  loading  spectrum  was  postulated  to 
be  a  random  variable.  The  sum  of  a  given  number  of  such  variables  represen¬ 
ted  the  total  crack  length  after  the  repetition  of  a  number  of  spectra. 
Failure  takes  place  when  for  the  first  lead  cycle  the  total  crack  length 
reaches  a  critical  value. 

We  assume  that  the  incremental  crack  extensions  per  cycle  are  indepen¬ 
dent  and  identically  distributed  variates.  If  X^  for  I  ■  l,...,n  are 

2 

the  increments,  each  with  mean  y  and  variance  o  ,  then  after  n  such 
cycles,  with  +...+  Xn  the  total  crack  length,  the  random  number 

N  of  such  increments  before  failure,  i.e.  before  reaching  a  critical  crack 
length  w  ,  must  satisfy 

'S  -ny 

n  >  a)-ny 

.  ''ll  a  /no. 

Hence  we  have  asymptotically,  by  the  central  limit  theorem, 

FN(n)  »  * 

where  is  the  standard  normal  distribution.  In  this  form,  this  distribu¬ 
tion  had  been  discovered  by  Freudenthal  and  Shinozuka  in  [  4  ] ,  before  its 


ny-u) 
Jn  0 


P[N  5  n]  =  P[S  >  w]  =  P 
n 
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publication  in  [  3  ] .  A  few  years  before  that  it  had  been  given  by  E.  Parzen 
[5].  Moreover,  T.  Von  Karman  had  found  it  in  the  twenties.  (One  is  led  to 
examine  the  collected  works  of  C.  F.  Gauss  in  this  regard.) 

The  difference  between  the  approach  in  [  3  ]  and  the  others  was  the  inter¬ 
pretation  and  parameterization.  Let 

0  ■  jj  ,  a  ■  o//\ Iw  . 

Then  the  distribution  can  be  written  in  the  form 


where  0  is  the  median  life  and  a  is  a  shape  parameter  and  where 

€(x)  *  Jx. - —  for  x  >  0  . 

/x 

Moreover,  estimation  procedures  have  been  derived  in  [  6  ]  for  this  dis¬ 
tribution.  This  distribution  is  related  to  the  log-normal  in  that  it  satisfies 
a  functional  relationship 

-  -?(i) 

which  is  a  relationship  between  a  variate  and  its  reciprocal.  We  now  ask: 

Is  it  possible  to  obtain,  by  other  simple  plausible  models  which  count  cycles 
until  failure,  some  useful  alternative  distributions  which  are  within  or 
similar  to  the  original  class? 

Let  us  first  suppose  that  the  critical  value  of  the  crack  length  at 

which  the  structure  fails  is  a  normal  random  variable  W  with  mean  u)  and 
2 

variance  p 
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Let  a  critical  value  w  be  given.  Then  the  number  N  of  cycles  to 
exceed  w  has,  by  the  preceding  argument,  the  distribution 


P[N  5  n|w]  =  1- 

Using  the  law  of  total  probability  we  find  the  distribution 


P[N  5  n]  -  f  0  SLSL  d 

J  w  1  p  ; 


=  * 


ny-ca 

v/na2+p2/ 


Let 


and  we  have 


where 


qu) 


5(t) 


for  t  >  0 


where  e  >  0  ,  small,  is  a  nuisance  parameter.  However,  note  that 


lim  £.(t)  =  »,  lim  5(t)  =  -  —  . 

t  "*■  s  t  -*■  0  i/c 

00 

This  value  determines  the  probability  that  the  structure  is  broken  (failed) 
at  time  zero.  (Such  an  event  has  been  known  to  occur  in  practice).  As  a 
consequence,  as  given  cannot  satisfy  the  reciprocal  relationship, 

to  wit,  C1(t)  =  -  Cx(“)  • 
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Let  us  also  examine  further  the  second  assumption  concerning  crack 

growth  made  in  [3];  namely,  given  the  total  crack  length  at  the  start  of 

the  loading  cycle  is  s  ,  the  incremental  extension  is  a  normal  variate 

2 

with  mean  p  +  6s  and  variance  a 

From  the  results  given  in  [  3 ]  we  know  that  the  random  number  of  cycles 
to  exceed  the  critical  length  oi  is 


p  -u 

FN(n)  -  *-§- 

n 


where 


I  (!«)*-!  I  ,  a*  -  o2fii±^i 

61  J  n  (1+6) 2-l 


If  6  >  0,  then  the  crack  accelerates  while  if  6  <  0  then  the  crack  de¬ 
celerates. 


Case  A:  6  >  0  . 


Let  ~  -  ln(l+6)  >  0.  Then 

p 


FN(n)  -  $ 


M1(en/6-l)-ai 


V*2e  _u,2 


where 


yl  6’  °1 


,  y2  -  — ,  id2 


'(1+5)  -1 


l+^-=l 

y 
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t 


Then 


Hence  we  have 


52<t) 


5(t) 


2t 

[C(eZt)]2 


for  t  >  0 


Again  we  note  that 


lim  £.(t)  ■  -  00 ,  lim  ^„(t)  ■  - 

t  +  0  L  t  +  00  /y" 


preventing  a  reciprocal  relationship. 
Case  B:  6  <  0  . 


Let  7T  “  -  ln(l+6)  >  0. 
P 


Then 


\  I 


P[N  <  n]  -  4> 


Pjd-e'"76)  -  0) 


where 


Let 


yx  -  ^  ,  a:  =  o//l-(l+6)2 


Pi 


Pj-w 


m2  o1  *  2  ax  * 


and  we  obtain 


P [N  <  n]  =  4> 


U!2en/,^-y2 
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which  is  exactly  the  same  model  as  in  the  preceding  case  except  that 
and  ^2  are  interchanged.  Thus  the  models  for  crack  acceleration  and 
deceleration  lie  in  the  same  family. 

These  arguments  demonstrate  that  it  may  be  naive  to  suppose  that  a 
more  realistic  theory  can  be  formulated  without  making  use  of  knowledge 
from  other  fields,  besides  mathematical  statistics,  in  the  construction  of 
distributional  models  for  fatigue. 


SECTION  III 


THE  RELATION  BETWEEN  GROWTH  RATE  AND  THE  DISTRIBUTION 
OF  TIME  UNTIL  CRITICAL  SIZE 

Consider  a  fatigue  crack  developing  in  a  given  structural  component 
due  to  the  repeated  impositions  of  a  loading  spectrum  through  its  service 
usage.  The  crack  growth  is  believed  to  be  divided  into  several  distinct 
phases,  each  one  of  which  is  governed  by  a  different  physical  mechanism. 
Neither  the  physical-metallurgical  basis  for  this  supposition  nor  the 
nature  of  the  mechanisms  involved  will  be  discussed  here.  In  this  regard, 
see  [  7  ]  and  [  8  ].  Instead,  an  appropriate  stochastic  model  will  be  pos¬ 
tulated  for  the  behavior  of  the  crack  in  each  phase,  each  of  which  consti¬ 
tutes  a  large  source  of  variation  in  the  observed  fatigue  life. 

The  first  phase  ends  when  the  crack  reaches  a  predetermined  detectable 
size.  This  is  called  the  initiation  phase.  We  label  this  minimum  detect¬ 
able  size  U)q  ,  which  is  a  function  of  both  the  material  and  the  inspection 
procedure.  The  time  of  initiation  will  vary  from  component  to  component 
even  though  they  are  nominally  identical  structures  subjected  to  the  same 
service  loads.  Let  T^  be  the  time  when  the  crack  reaches  the  specified 
initiation  size  .  We  assume  that  T^  is  a  non-negative  random  variable, 
across  the  population  of  components,  with  a  distribution  which  is 

unspecified  at  present:  T^  ~  .  We  do  not  exclude  the  possibility  that 

Tq  has  a  positive  probability  of  being  zero.  This  would  correspond  with 
the  situation  when  the  specimen  was  initially  flawed. 

During  the  second  phase  of  crack  growth  the  crack  grows  at  a  rate 
which  is  virtually  constant  for  each  component  and  loading  spectrum  examined 
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but  the  rate  will  be  somewhat  different  for  different  specimens  even  if  the 
same  loading  spectrum  is  imposed.  Let  s(t)  denote  the  observed  crack 
length  at  any  time  t  ?  0  .  Given  that  the  crack  was  initiated  at  time  t^  , 
our  assumption  is  that  for  a  time  thereafter  the  growth  rate  is  constant, 
namely 

s' (t)  =  u  for  t  >  tg  . 

Here  the  slope  u  is  determined  by  the  material,  the  work  cycle  and  the 
specimen  itself.  It  now  follows  that  the  observed  crack  length  would  be 
given  by 

s ( t )  *  u(t-tQ)  +  (*)  for  t  >  t0  . 

The  constant  growth  rate  of  the  crack  observed  during  the  second  phase  of 
behavior  varies  from  one  metallic  specimen  to  another  and  thus  the  stochastic 
nature  of  the  crack  follows.  The  population  of  crack  length  of  these  metallic 
specimens  can  be  described  by  a  stochastic  process,  given  the  crack  initiation 
occurred  at  time  t^  ;  to  wit , 

S(t)  *  u(t-tg)  +  Wg  for  tQ  5  t  2  tx  (3) 

where  U  >  0  is  a  random  variable  with  a  distribution  yet  to  be  specified 
across  the  population  of  specimens. 

Because  of  the  multitudinous  influences  governing  each  incremental 
crack  extension  the  crack  length  at  a  given  time  can  be  regarded  as  a  Gaussian 
Stochastic  Process;  namely, 

S(t)  ~  31  (y t ,  a2)  for  each  t  >  tQ.  (4) 
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Moreover,  for  t  >  tg,  we  have 


°0  ^ 


and  p  “  ml (t-t0>  +  W0  * 


From  equation  (4)  it  follows  that  the  probability  that  at  cime  t  >  tg  the 
crack  does  not  exceed  a  specified  length  to  >  Wq  is 


P[S(t)  >  u|T0  =  tQ] 


for 


t  >  t 


0  ' 


(5) 


Here  the  constant  m^  is  the  expected  rate  of  crack  growth  as  determined 
by  the  loading  spectrum  and  the  material.  The  time  scale  is  chosen  so  that 
each  application  of  the  spectrum  occurs  in  unit  time.  We  do  not  regard  it 
as  an  unknown  parameter  to  be  estimated  statistically.  It  should  be  calcu¬ 
lated  from  the  loading  spectrum  and  the  material  by  the  use  of  a  formula 
such  as  the  one  given  in  (  9  ] . 

Let  Tj  be  the  time  at  which  the  crack  first  reaches  the  critical 
propagation  length  oo^  and  stable,  linear  growth  ends  and  the  third  phase 
of  growth  begins.  During  this  phase  a  differential  equation  following  the 
laws  of  fracture  mechanics  is  presumed  to  describe  the  behavior.  The 
critical  length  is  again  a  constant  calculated  from  the  work  load  and 

the  material. 


Thus 

Define 


inf {t:  S(t)  > 


^0 

m. 


/m^ (Wj-Wq) 


Then  from  equation  (5)  it  follows  that  for  t  >  tg 
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where 


P‘T1  ‘  ‘lTo  ■  ‘o1  ■ 


5(x)  - 

Thus  we  see  that  the  distribution  of  is  ^-normal  with  characteristic 

value  3^  and  shape  parameter  . 

But  by  suitable  rearrangement 

PITj  5  1 1  Tq  -  tQ]  *=  P[U  ?  (u>1-a»0)/(t-t0)3 

and  it  follows,  since  £(l/x)  =  ~C(x)  ,  that 

P[U  5  u]  H  £(u/m^)j  u  >  0  • 

Thus  the  distribution  of  the  growth  rate  (slope)  U  of  the  crack  during  its 
stable  growth  is  also  ^-normal  with  characteristic  value  m^  and  shape 
parameter  .  This  reciprocal  property  of  the  ^-normal  (Birnbaum-Saunders) 
distribution  has  been  discussed  in  [  6  ]  and  []_()]• 

We  summarize  this  discussion  with  the  following 
Theorem  1:  Following  crack  initiation  of  size  uig  at  given  time  tg  , 
if  the  crack  length  S(t)  for  time  t  >  tg  is  a  Gaussian 
Stochastic  Process 

S(t)  ~U(ut,o^) 


P[s(t)  ;  Ul|T0  -  t0l 


t 


“l  V  8! 


/x -  for  x  >  0  . 

& 


=  °0^ 


m^t-tg)  +  Wg 


for  t  >  t. 


* 


where 


and  each  crack  progresses  at  a  uniform  rate  until  a  critical  propagation 
size  is  reached  so  that 

S(t)  -  U < t-- tQ )  +  u>Q  for  t  >  tQ 

then  the  random  rate  of  progress  U  is  a  ^-normal  with  characteristic 
value  ra^  and  shape  parameter  * 

U  ~  ^(a^n^)  • 

Moreover  the  time  ,  following  initiation,  until  the  critical  size  u)^ 
is  reached  is  also  ^-normal  with  characteristic  value  3^  and  shape 
parameter  , 

T1  ~  9lC(a1,61) 

ao  „  vwo 

-  ,  K  *  ~  . 

"l^i  1 


where 
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SECTION  IV 


THU  DISTRIBUTION  OF  TIME  UNTIL  FATIGUE  FAILURE 


During  the  third  phase  of  fatigue  crack  growth  the  crack  is  assumed  to 
be  extended  as  a  function  of  a  propagation  factor.  This  propagation  factor 
is  a  function  of  time  which  has  been  taken  as  the  solution  of  a  differential 
equation.  This  solution  was  derived  from  an  assumption  taken  from  fracture 
mechanics:  namely  that  the  rate  of  propagation  of  the  crack  is  a  linear 
function  of  the  square  root  of  the  crack  length,  see  [11]. 

For  given  values  of  a,b  >  0  ,  the  propagation  factor  is  of  the  form 


p(t:  a,b) 


eat  if  b  -  0 

(l+abt)1/b  if  b  *  0  . 


If  we  suppose  that  b  is  fixed  for  a  particular  geometry  and  material  then 
the  controlling  constant  a  will  have  a  distribution  across  specimens  drawn 
at  random  from  the  population  of  such  components. 

It  is  then  clear  that  the  time  spent  in  the  third  phase,  say  T^  ,  is 
given  by 

7  £n(a3_-a)0)  if  b  -  0 

A  J  Z 


T 


3 


u 


(u>3-u>2) 

b“ 


if  b  *  0  . 


Thus  if  we  assume  A  is  ^-normal  then  it  would  follow  that  T^  is  ^-normal 
as  well. 

Hence  the  time  until  fatigue  can  be  represented  as  being  the  sum  of 
three  random  quantities:  the  time  until  crack  initiation,  the  time  between 
initiation  and  critical  size  and  the  time  between  critical  size  and  ultimate 
size. 
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Crack 

length 


Ultimate 

size 


Critical 

size 


Initial 

size 


The  relative  proportion  of  time  which  a  crack  spends  in  the  three  various 
phases  can  not  be  determined  in  practice.  However,  there  are  instances  where 
the  crack  initiation  occurs  very  late  and  others  where  it  was  suspected  to 
have  been  initially  flawed.  It  is  also  thought  that  the  time  spent  in  the 
third  phase  is  relatively  short  compared  with  the  other  two. 

Each  of  these  distinct  phases  has  been  studied  in  the  laboratory  not 
only  for  different  materials  under  different  environmental  conditions,  but 
different  geometrical  configurations.  Despite  this  extensive  effort  there  is 
still  so  much  variability  in  material,  workmanship  or  service  exposure 
that  no  calculation  of  fatigue  reliability  based  upon  a  synthesis  of  these 
theories  is  universally  accepted  prior  to  service  with  the  same  assurance 
that  the  final  fatigue  test  carries. 

Now  we  come  to  a  critical  point  in  the  development.  Suppose  that  we 
are  interested  in  the  convolution  of  non-negative  random  variables  for 
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i  ■  l,2,...,k  which  are  independent  but  not  necessarily  identically  distri¬ 
buted  by  one  of  the  £ -normal  distributions.  In  particular  we  think  of 
T1  +  T2  +  T3  where 

is  the  time  from  the  introduction  of  service  until  crack  inititation; 

T2  is  the  time  between  crack  initiation  and  critical  crack  length; 

is  the  time  following  critical  crack  size  until  crack  growth  terminates. 
Here  crack  growth  termination  may  be  at  component  fracture  or  when  the  crack 
reaches  a  "crack  stopper"  which  is  located  in  certain  components. 

We  ask,  what  is  the  distribution  of  the  total  time  until  crack  termina¬ 
tion? 

Here  and  subsequently,  iff  means  if  and  only  if. 

Recall  that  T  is  ^-normal  iff  for  some  d,B  >  0 

Z  *  £  5(T/6)  ~  31(0,1)  , 

where  £  is  a  monotone  increasing  map  of  (0,°°)  onto  (-00,00) 
such  that 

£(x)  =  -  £(l/x)  for  x  >  0  . 

This  means  that  tp  -  £  ^  must  satisfy  the  relation 

«K-z)  =  for  -oo  <  z  <  00  and  ^(0)  =  1  .  (6) 

We  contend  the  following:  If  are  sufficiently  small  and  are 

standard  normal  variates  then 

£  Bi^i(aizi)  =  (I  8i)^(£  a^)  (?) 

for  some  function  ip  which  satisfies  (6)  as  .  ^  . 
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Note  that  for  a sufficiently  small  we  can  certainly  make  the  approxi¬ 
mation  given  in  (2)  using  the  first  few  terms  of  the  Taylor's  series  approxi¬ 
mation.  We  now  note  that 


Z  ViVi’  ■  I 


V-Vi* 


i  <r  Bj)2  [I  e1«’1<-a1zi)l"1 


(8) 


The  last  approximation,  actually  an  inequality,  follows  from  the  relationship 
between  arithmetic  and  harmonic  means  and  is  exact  when  all  but  one  of  the 
B^'s  are  zero.  Thus  by  applying  the  same  argument  as  before,  namely  that 
(7)  must  hold  with  z^  replaced  by  -z^  ,  and  then  assuming  (8)  is  exact 
we  have 


and  hence 


Z  6iipi(a1z1) 


E  6i 

iK-Zaizi) 


ij;(Zaizi) 


1 

iK-Za^) 


We  now  claim  that  if  T^  are  ^-normal,  then  as  a  mathematical  approxi¬ 
mation 

T  ■  T^  +  T^  +  T^  is  also  ^-normal 
but  the  £  is  different  from  for  i  =  1,2,3.  Moreover, 


we  claim  that  as  an  empiiical  fact  this  approximation  is  supported  by  the 
evidence  which  has  been  accumulated  on  the  total  time  of  service  from  intro¬ 
duction  until  fatigue  fracture.  Moreover  we  feel  that  ^-normal  variates, 
which  are  close  to  log-normal  in  construction,  are  sufficiently  descriptive 
of  the  fatigue  phenomenon  that  their  investigation  is  of  practical  importance. 
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Thus  we  proceed  with  the  subsequent  examination  of  the  problems  of  estimation 
which  arise  in  this  connection. 

The  final  fatigue  test,  which  can  be  completed  on  only  a  few  specimens, 
because  of  the  time  and  expense,  confounds  the  effect  of  each  separate  phase 
of  the  fatigue  crack  growth  within  the  observation  of  the  total  life. 

Since  the  number  of  complete  specimens  which  can  be  used  for  observations 
of  the  fatigue  life  is  so  small  an  extremely  efficient  theory  of  estimation 
is  needed  which  is  specifically  tailored  to  the  data  which  can  be  obtained, 
both  in  quantity  and  form. 
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SECTION  V 


THE  CLASS  OF  DISTRIBUTIONS 


Let  E  be  a  class  of  real  valued  functions,  concave  monotone  increasing 
and  mapping  the  positive  real  line  onto  the  real  line  such  that  if  £cE 
then  it  satisfies  the  functional  equation 

$(t)  -  -£(l/t)  for  t  >  u.  (9) 


If  T  is  a  non-negative  random  variable  such  that  there  exist  con¬ 
stants  a,p  >  0  and  £eS  for  which 


£  C(T/8)~  91(0,1)  (10) 

we  say  that  T  is  Z-nomal;  see  (10j.  This  nomenclature  is  adopted  by 
analogy  from  that  of  the  log-normal  distribution. 

Thus  if  T  is  ^-normal  with  parameters  a,B  >  0  and  has  distribution 
F  then 


F(t) 


7  «t/S) 

0 


for  t  >  0 


where  $  is  the  standard  normal  distribution  given  by 

*(*)  -  —  |  e"t2/2  dt. 

/27  j-°° 


Virtually  all  the  collections  of  fatigue  life  data  have  been  recorded 
in  terms  of  the  logarithm  of  the  number  of  cycles  until  failure.  Moreover 
these  data  have  been  plotted  on  log-normal  probability  graph  paper.  To 
facilitate  the  comparison  of  these  new  families  of  life  distributions  with 
plotted  data  currently  available,  we  shall  consider  the  log-life  variate 
X  =  inT  and  we  take  its  distribution  G  to  be  of  the  form 
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•jnnij*.  :  'mm  WW  (MM 


where 


G(x)  -  « 


-  w(x-y) 


-OD  <  X  <  00 


p  «  in{5,  to(x)  -  £(ex)  for  all  x. 


(ID 


We  now  consider  the  class  of  functions  ft  defined  as  follows:  For 
a  real  valued  function  u,  mapping  the  real  line  onto  Itself,  which  is  twice 
differentiable  we  say; 

weft  iff  a i*  ^  0,  to'  >_  to"  and  to  is  odd.  (12) 

Let  us  examine  the  closure  properties  of  ft. 

We  state  without  proof: 

Theorem  2:  If  t^.o^  e  ft  then  to^  +  <o2  e  Moreover  if  to  e  ft  and 
a>0,  0  <  b  <  1  then  both  to^u^  e  R  where,  for  -»  <  x  <  °°  , 

to^x)  *  ato(x)  and  w^x)  ■  <o(bx). 

As  a  consequence,  for  any  given  to  e  ft  we  can  generate  a  three- 
dimensional  parametric  family  of  distributions  of  the  log-life,  namely, 

G (x : a , y , p )  -  sjjio^J^J  -«  <  x  <  00  (13) 

where 


i 

l 


a>0,  y  >  0,  — 00  <  p  <  °°  , 

In  this  notation  p  is  the  location  parameter,  y  is  the  scale 
parameter  and  a  is  the  flexure  parameter.  The  flexure  and  scale  together 
control  the  shape. 
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Examples  of  such  functions  <u  which  correspond  to  known  parametric 


families  are: 


u>^(x)  ■  sinh  x  -«  <  x  <  «®  (13) 

and 

uijCx)  -  $  ^[F(x)J  -oo  <  x  <  00 

where 

F(x)  »  —  +  (sgn  x)  exp{-(£n  2)e  ^}. 

Here  u>^  corresponds  to  the  distributions  of  Birnbaum- Saunders  in  [  6  ] 
and  ujj  corresponds  to  the  symmetric  Weibull  distributions  of  [  12] . 

There  are  several  ways  of  constructing  functions  in  ft.  One  is  as 
follows:  Given  any  real  valued  function  f,  mapping  the  real  line  onto 
itself,  such  that 

f’  >_  f"  >  0  define  w(x)  ■  f  (x)  -  f(-x). 

A  second  method  follows  from  Theorem  2,  For  example,  if  we  have  constants 
a^  >  0,  1  >  >  0  for  i  *  l,...,n  then  weft  where 

n 

u(x)  =  2a^sinh(b^x)  <  x  <  »  , 

1 

Now  we  state  a  result  on  closure  under  inverse  transformations  in 

2  -1 

Theorem  3:  If  ^eft  and  (w1)  >_  -u"  then  w  eft  . 

Proof:  Clearly 

(w  ^)  ’  =  l/u'(w  ^)  >_  0  since  to'  >_  0, 
and  to  odd  implies  to  ^  is  odd.  We  must  only  check  that  (a)')"  <  (to  ^)  * 


22 


By  definition 


[u’to”1)!2  [u'Oa”1)]3 

But 

-  “"fa'1)  <  1 

[u'Cu"1)]3  “■•(■T1) 

is  equivalent  with 

-u"(w“1)  <  (UJ-1)]2  iff  -0)"  <  (U)')2  .  || 

Thus  there  are  inverses  of  certain  elements  of  ft  which  are  themselves 
Included  in  ft.  Using  the  notation  of  (13)  we  have 

Theorem  4:  ui^eft  where  ui^Cx)  ■  sinh  x,  and  we  define 

w3(x)  -  w^Cx)  «  £n(x  +  /xZ  +  1  )  .  (14) 

Proof:  By  definition 

«  sinh,  ■  cosh,  ui^  ■  w^. 

We  must,  by  theorem  2,  check  that 

2 

-sinh  x  <_  (cosh  x)  for  all  x  ,  -<»<  x  <  <»  , 

If  we  let  x  be  replaced  by  -x,  this  inequality  is  equivalent  with 

2 

-sinh(-x)  <_  (cosh(-x)]  . 

But  by  known  properties  of  the  hyperbolic  functions  this  is  equivalent 
with 

2 

sinh  x  <  (cosh  x)  -«  <  x  <  «  . 
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From  the  definitions  we  must  show  that 


(ex-e“x)/2  <_  (e2x+2+e~2x)/4  , 


which  is  equivalent  with 


-  2x  _  x  ,  o  .  0  -x  .  -2x  ,  x  ,,2  ,  ...  -x,  2 

0  <_  e  -2e+2  +  2e  +e  -  (e  -1)  +  (1+e  ) 


Let  us  illustrate  the  typical  behavior  of  one  such  family  and  its 
reciprocal.  Consider  the  two  three-parameter  families  induced  by  (13) 
for  u  >  0,  y  >  0,  -  oo  <  y  <  oo  which  are  generated  by  and  uy 


“(x)  *  l  sinh  (^) 


_  1  1  o_  r  x-u  +  /(x-u)^+  Y7  1 


Figure  2.  Specific  transformations;  x  =  y+z  plotted  against  u>(z)  and  $ 


SECTION  VI 


25 


3L  3L  3L 

Now  the  joint  solution  of  —  »  0,  —  ■  0,  ■  0  will  yield  the  maximum 

likelihood  estimators  a,Y,y.  We  thus  seek  the  simultaneous  solution  to 
the  following  equations  in  the  variables  (a,Y»u)s 

i  ?  /xru\  ./V\  ?  "I— ) 

7  |  I  •  <19 

n  x  -yv  /X.-Wv  ,x.-Pv  ,  n  ,x.-y  x  «"(■— — } 

ii  y I  y  -L  — 1JL.Z  ,  i  ,  (20 

a2  n  j-  \  y  /  \  y  f  \  y  '  n  ^  \  y  /  /xi"M\ 

“'(■7-7 


o2  - 1 !  *2(^) . 

nr  \  Y  / 


As  a  notatlonal  convenience,  let  us  define  as  implicit  functions  of 
y  and  y  the  values 


f±  y 


i  ■  l,...,n 


and  introduce,  for  any  function  f,  the  averaging  operator 


<f(’rl»  ‘  n  I  Hvl>  ■ 


We  define  as  functions  of  (u,y): 


_  <P(y1» 

H  -  - <R(y1)>  , 

<« (y±)> 


<yiP(y1)> 

<t»2(yi)> 


-  <yiR(yi»  -  i  » 


where  we  have  the  ratio  and  product  functions  defined  when  both  R  and  P  are 


odd,  by 


R  ■  ,  P  »  u)  u) *  . 
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Let  us  suppose  throughout  the  following  discussion  that  we  have  several 
samples  from  different  populations.  Suppose  we  have  m  such  samples 


x.  . »  • • • »x  a  for  J  “  l»...»m. 
nj'J 

where  n^  observations  are  taken  from  the  population  with  parameters 
for  j  ■  l,...,m.  Let  us  compute  the  two  sample  statistics 

"j  "j 

x.  *  —•  y  x  ,  B2  -  ~r  V  (x  -  X  )2 
J  nj  iJ  J  nj-1  itl  iJ  j 


X. .  ~  I(°t4  »Vj  tU.)  for  i  -  1, . . . ,n. ; 

J  j  J  J 


ai  \  yi  ' 


1  "  1 . nj;  j  = 


are  independent  identically  distributed  standard  normal  variates.  Since 


a)  is  odd  we  see 


y  3^  -  E  u~\a 


JV  - 0 


and  hence 


E  Xij  =  Mj 


and  thus  Xj  is  an  unbiased  estimate  of  .  But  note  also  that 

r  i  ii  /X.  .-y.  \  2 


e[“'1<vu,]2-e(-V)  • 
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If  we  define 


°)  ’  var(V 

we  see 

2  2  _  -1,  7  .2 

“i  ■  E "  (“jzij) 

2 

and  we  know  that  an  unbiased  estimate  of  is  the  sample  standard  devi¬ 
ation  Sj. 

Thus  we  see  that  this  family  of  distributions  separates  the  variance 
into  the  product  of  two  factors.  The  first  factor  is  the  scale  parameter 
while  the  second  factor,  say  B(a),  is  determined  by  the  flexure  parameter 
a  namely 

B(o)  -  E  u)_1(aZ)2  -.L0u"j( y>2d*(^)  • 

Thus  the  equation,  valid  for  large  samples,  of 

s3  ■  A  b(v 


forces  us  to  utilize  an  independent  method  for  estimating  a 
We  now  prove  the  important 


y 


Theorem  5:  The  maximum  likelihood  estimator  of  ,  defined  by 


nj  i-1  V  YJ  > 


-i  *5 


has  a  distribution  independent  of  and  y 


(23) 


Proof:  Let 


'tj 


!u2i 


for  i  ■  1, . . . ,n 


be  standardized  variables  independent  of  p^  and  y  and 
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define 


U 


J 


Yj 


V  ■  Y  /y 

J  YYJ 


as  the  normalized  maximum  likelihood  estimators.  We  now  write 

Y  -U. 


T 

1J 


J 


tm  -  ^ 


(24) 


and  we  note  that  for  each  j  the  m.l.e.'s,  namely  and  y^ , 

satisfy  the  two  equations  which  we  can  write  as 


<P(T1;j)>  ■<R(Tlj)X»(Tij)> 


<Tijp(Tij)> 


<TijR(Tlj)>+  1 


<«^(T1J)>  , 


making  use  of  the  product  and  ratio  functions  P  and  R  as  defined 
in  (22)  .  The  averaging  operator  <•>  Is  taken  only  over  the 
subscript  1.  Since  the  for  1  ■  l,...,n^  have  distributions 

Independent  of  vij  and  y ^  ,  it  follows  a  fortiori  that  so  must  a^. 


Presume  that  we  have  computed  the  maximum  likelihood  estimates  previously 
deflnevl  for  m  different  groups  of  data  as 


A  A 


V  Y  UJ  for  i  * 1 . "• 


Under  the  null  hypotheses  that  all  the  aj'8  are  the  8ame» 

H0:  “l  ■  °2  ■  •••  ■“»  • 


we  might  take  a  grand  average 

-2 
a 


m 


m 


I  ».  K/  i  ». 

j-1  J  j-1  J 


as  the  combined  estimate. 
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We  know  from  general  statistical  theory  on  the  consistency  of  the 

*2  2 

maximum  likelihood  estimators  that  -*>  in  probability  as  Oj  -*■  ». 
Unfortunately  we  have  rather  small  sample  sizes  n^  for  j  ■  l,...,m  but 
the  number  of  groups  m  may  be  quite  large.  Ihus  if  it  were  possible  we 
would  determine  bias  factors  b(n),  as  a  function  of  the  sample  size 
n  “  1,2,...,  such  that 

E  -  .J  b^, 

from  which  we  could  define  unbiased  estimates  and  have  consistency  as 
m  as  well.  However  from  (3.1)  we  have 


e  yE  cj2^) 

where  T^  is  defined  in  (24).  Thus  one  cannot  guarantee  that  can 

be  factored  out  leaving  b  a  function  only  of  the  sample  sizfl.  We  con¬ 
jecture  that,  in  fact,  it  is  the  case  that  b  will  also  be  a  function 
of  a. 

We  now  turn  to  the  problem  of  estimating  the  common  a  and  the  related 
problem  of  testing  the  hypotheses  Hq  from  the  point  of  view  of  general 
statistical  theory.  Let  us  suppose  that  Hq  is  true  so  that 


ij 


I(a,Yj iMj) 


for  i  “j  • 


The  joint  likelihood  is  the  sum  of  likelihoods  for  each  of  the  m  sets 
of  data,  each  one  of  which  is  similar  to  the  one  before,  namely 


m 


L  -  I  L 
J-l 


j 


where  by  comparison  with  (  15  )  we  see 
n, 


“J 

£ 


-in(ay^)  - 


+  ii m 


*4  4'^ 


(0 


2a 
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By  comparison  with  (16)  wa  have  for  j  ■  l,...,m 

it-  .  y  _1_  pto)  .  i  R  to) 

8“j  1^1  .‘r,  \  TJ  )  'i  \  TJ  / 

where  we  have  made  use  of  P  and  R  as  defined  In  (22)  .  Similarly, 
we  have  from  (  17 )  for  j  ■  1  9  •  •  •  |  ID 


and  from  Eq.  (18) 


SL  "  n>  r  ^  1  2/Xlt'Ul\l 

- •  *  ^  "7  ■  nnj • 


The  method  of  obtaining  the  maximum  likelihood  estimators  then  is  as  follows: 
We  must  solve  for  the  2m  +  1  variables 


«.Yi,...Yn.Mi,...,Mm 


in  the  2m  +  1  equations  which  we  can  write,  letting  y^ 
"j  nj 

if  l  p(yn)--r  1  r(Yh>  j-i . ® 

a2  nj  1-1  «  nj  iti  « 


7 1  i,  yii  "v  - 1  Jx  ^  R(j,u)  * 1  j  ■ 1 . ”■ 


n, 

®  j  O  ,  « 

1  I  ni 

j-1  1-1  3  j-1  3 
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The 


*  *  A  A  A 

Let  us  call  the  simultaneous  solution  a  ,y  ,  ...,y  . . .  . 

mi  m 

actual  method  of  solution  we  defer  until  later. 

Let 

6  "  (“l . VY1 . VW1 . Vm)m 

Define  the  parameter  space  and  a  subset  Qq  as  follows: 

0  e  Iff  >0,  >  0,  -oo  <  <  *  for  f  .  l,,.,,m 

0  e  iff  0  e  !1,  and  a,  ■  a~  *  . . .  ■  a  . 

0  1  12  m 

These  two  sets  correspond  to  the  two  hypotheses. 


Hq:  6  e  Aq,  0  e  . 

The  likelihood  ratio  is 

8up  n  n 

6eii0  ^-x 

"  sup  n  n 

6eQi  j-i  i-i^VYYV 

We  know  by  statistical  theory  that 

n. 

m  J  it  it  it  *  *  * 

-2£nX  »  -2  V  £  [«,n  gU^ia  ,Yj*Pj)  -  &n  g(x^  :a^  ,Yj  ,Wj )  ] 


will  have,  asymptotically,  a  chi-square  distribution  with  3m-(2m+l)  ■  m-1 

degrees  of  freedom.  Thus  the  critical  region  will  be  the  upper  tail  of 
2 

the  x  (m-1)  distribution. 

In  order  to  compute  this  test  statistic,  we  use  (  15  )  and  set 

A 


yu 
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r 


to  write 


»  "J 


-2£nX  -  -2  S  S 
j-1  i-1 


*  *  1  2  *  * 

-  £n  a  y, - t-j  “  (y4i>  +  £n  “'(y,.,) 

3  2  (a  V  J 


+  in  a  y.  +  — 7 — r  w2(y  )  -  in  u)'(y  ) 

j  J  of*  \z  lJ  *0 


2<V' 

This  is  the  statistic  that  oust  be  computed  to  perform  the  test. 
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SECTION  Vil 


THE  BEHAVIOR  OF  THE  MAXIMUM  LIKELIHOOD  FUNCTION 
IN  CERTAIN  CASES 


We  now  consider  the  behavior  of  the  likelihood  function  when  co  is 
"s-shaped"  in  its  behavior.  The  typical  shape  is  that  of  w*  as  graphed 
in  Figure  2.  Formally  we  say  that  <*>  is  s-shaped  if 

(i)  1  >  >  u)'(x)  >  0  ,  03 *  (0)  =  1 

(ii)  u)' (x)  decreases  as  |x|  increases. 


Let  us  consider  the  joint  log-likelihood  equation  ii  Lving  only  the 
two  parameters  u,Y  ,  having  previously  eliminated  a  by  substitution. 
We  obtain 


LO-.Y)  •  -  -Vj  £  “2(-^— )  ♦  ±£  in 

2"fa>  i-1  i-1 


(25) 


where  a  is  a  function  of  y,y  given  by 


to2  •  iZ'AT)  • 


i-1 


We  find  after  substitution  and  simplification  that 


L(u,y)  =  -  y  In 

11  v  -II  2 

i  y[Y  J  i  y\| 

n  L  |Y  a,l  Y  /I 

1 

-  ?  +  ;  2>  “'(-T")  •  (26) 
1 
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uk  ki  .in  — «wui m  K.iM«Muwit>wi?  -witiBwWWj  WWSdlW’iBI'I#  ife? 


Let  us  set 

n 

y  *  x  -y  for  i  ■  l,...,n  ,  a  ■  max  |y  |  . 

1  1=1  1 

Now  from  property  (i)  we  have  for  i  =  1 . n 


[yw(yi/y)]2  >  y^  [u> '  (y±/y)  ] 2  >  y2[to'  (a/y)  ]2  . 


Upon  substitution  we  find 


L(li,y)  +  -  in 


1  z 

n  i—* 


+  -  <  -  £n  to'  (a/y)  . 


On  the  other  hand  we  have,  by  the  bound  in  (i) 


ty  t)(yi/y)]2  =  y\ 


w(y1/y) 


5  y^  for  i  =  1, . . .  ,n 


and  thus 


L(y,y) 


£n 


1  V  2 
n  2-  *i 


“  7  +  i  Z  £n 

i=l 


From  (ii)  we  know  that 


u>'  <y±/y)  >  u'(a/y)  for  i  =  1 . n 


(27) 


and  thus  from  (27)  we  have 


L(li,y)  +  £n 


+  j  >  £n  to'  (a/y) 


We  can  now  state 

Theorem  6:  The  log-likelihood,  as  a  function  of  y  >  0  ,  is  bounded  by 
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+  -|  5  -  Jin  oi'  (a/y) 


Jin  oj *  (a/y) 


L(y,y)  +  -  Jin 


L1 

n 

i-1 


(xry)‘ 


where 


n 

a  *  max  |  x .  — VJ  | 
i=l  1 


It  follows,  of  course,  that  in  the  limit  as  y  -*■  00  the  log- likelihood 
approaches  a  function  of  y  which  has  a  maximum  at  the  value  y  =  x  .  More¬ 
over,  certain  of  those  go's  chosen  for  numerical  examination  show  that  for 
fixed  y  ,  the  likelihood  is  an  increasing  function  of  y  >  0  and  moreover 
that  the  limit  is  reached  very  quickly.  To  the  postulated  behavior  of  s- 
shaped  go's  in  Q  we  add  the  assumption 


(iii)  Z  (x)  =  X  U decreases  as  I x I  increases, 

oi  0)(x)  1  1 

For  the  remainder  of  this  discussion  we  fix  y  ,  omit  its  mention  in 
equation  (26),  and  write  for  the  likelihood  simply 


L(y) 


where 


hj(y) 


Jin  oi'  (y±/y ) 


and 


y^  =  x^-  y  for  i  =  1 , . . .  ,n 

hi(y)  =  y  w(|yi|/y)  . 


This  equality  follows  since  oi  is  odd  and  positive  for  positive 
argument  and  therefore  xoi(x)  and  oi'(x)  are  even  and  positive  for  all 
real  x  . 
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li — it~r — - ^yvtMTv 


By  adding  and  subtracting  the  quantity 
combining  terms,  we  obtain,  in  new  notation. 


1 

n 


X  $,n  h  (y) 
1  1 


and  then 


l(y)  -  c(y)  +  d(y)  -  y  -  i  X  MyJ 

1 


(28) 


where  the  functions  C  and  D  are  defined  respectively  for  y  >  0  by 

n 

c(y)  -  £  X  £n  ^lyj/v)  (29) 

and 


D(Y) 


in  hi(Y) 


(30) 


We  study  their  behavior  in 

Theorem  7:  The  function  D  is  always  negative  and  decreasing,  while  C 
is  negative  and  increasing. 

Proof:  To  see  that  D  is  negative  we  recall  Jensen's  Inequality,  that  for 
any  random  variable  X  ,  if  0  is  convex  then 

4>(EX)  5  E  <f)(X)  . 


Since  -  —  in  x  is  convex  (and  decreasing)  we  see 


Jin  h , 


1 

from  which  we  conclude  that  D  5  0  . 


1 

To  see  that  D  decreases  we  show 
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This  is  equivalent  to 


where 


and 


aibi 


a±  -  hJ(Y) 

YhJ(Y) 
bi  “  h1(Y) 


[Y  ^(lyJ/Y)]2  >  0 

l  -  chJ(|yil/Y)  >  o  . 


(31) 


But  the  inequality  (31)  is  true  provided  that 


(ai~aj )  ^bi"*bj )  -  0  for  a11  • 


(32) 


(See  Hardy,  Littlewood  and  Polya,  [13]  p.  43) 


.  Firstly  since  to  (x)  is  an  even  function 
which  increases  as  |x|  increases  we  have  for  fixed  y 


Suppose  |yj  <  ]y 


=  y2  w2(|yi|/Y)  <  Y^  w2(|yj|/y)  = 


Secondly  since  £(x)  is  also  an  even  function  but  one  which  decreases  as 
]  x |  increases  we  see  from  (iii)  and  our  supposition  that 

bi  -  1  “  ^(lyJ/Y)  <  1  "  ylYjl/Y)  =  bj  . 

Thus  we  see  the  sequences  are  similarly  ordered  and  (32)  is  true.  Hence 
D'  is  negative  and  therefore  D  is  decreasing. 

To  see  that  C  is  negative  we  note  from  (iii)  that  1  *  £^(0)  >  C^(x) 
and  since  £  decreases,  it  follows  that  Jin  ( I y ^ I /y)  increases  as  a 
function  of  y  .  |  | 
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We  now  show  that  D  is  virtually  constant. 


Theorem  8:  Let 


Vy)  “  lim  for  y  '  0 

X  ■>  00 


D(Of) 


2S  E 

i-1  j-1 


Moreover 


n  n 

“<")  ■  £E  ‘“lyj  ■  t  ln  jE  i^il2 

1  j-1 


D(°°)  <  D(0+)  5  0  . 


Proof:  Note  that  h^y)  *  |  |  <t>  (  |  y  ^  |  /y)  i  =  l,...,n  and  since 

<l>(0+)  *1  we  have  lim  h  (y)  *=  |  y  |  .  From  this  fact  and  the  defini- 

Y  -►  oo  i 

tion  given  in  (30)  we  see  that  D(Q0)  is  as  given.  An  expression  equivalent 


to  (30)  is 


n  f  n  ,  2 ,  x 

riln  iI4^- 

'"w 


Now  by  definition: 


lim  [h  (y)/h. (y)] 

Y  -*■  0  J 


<*)(|yj/yilx) 


X  -*■  00 


njiy./yil) 


and  hence  the  result  follows. 
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We  now  show  that  the  range  of  C  is  quite  extensive  . 


Theorem  9: 

In  all  cases  C(°°) 

■  0  ,  while 

in  6  if  6  >  0 

C(0+)  =  ■ 

.  -  00  if  6-0 

where 

6  =  lim 

6  (x)  0  5  {  < 

to 

1  • 

x  -*■  “ 

We  now  show 

the  relationship  of  the  S-shapedness  of  a  function  and 

the 

behavior  of 

its  inverse: 

Theorem  10: 

to  is  S-shaped  i.f 

it  satisfies  properties  (i),  (ii) 

and 

(iii)  iff 

its  inverse  <p  =  to 

1  satisfies  (i) * , 

(ii)’,  (iii)’  where 

(i)' 

!  5  iUii.  <  u'(x>  , 

(ii)' 

<K(x) 

increases  as  |x| 

increases  , 

(iii)' 

x^' (x) 
'i'(x) 

increases  as  |x| 

increases. 

Proof:  First  note  that  u)  =  ip  and  to'  (\p)  =  ,  then  let 

x  =  ip(y)  ,  which  is  an  order  preserving  transformation,  so  that 

to(x)  =  <4iHy)  ]  =  y  'P'  (y) 

x  u)'  (x)  <Ky)  to ’  [4^ (y )  ]  (y ) 

Thus  we  see  the  claim  is  true  for  (iii)  '  .  The  others  are  checked  in  a 
similar  way. 


AO 


Consider  the  boundary  case  where 


U)(x)  ”  |x|asgn  X  -  00  <  x  <  00 

for  some  a  >  0  .  Now  a)  satisfies  property  (iii)  iff  a  <  1  .  But  it 
follows  that 


a) 


Consequently  from  the  form  of  the  likelihood  given  in  (28)  we  see  that  C 
is  constant.  Moreover  we  check  easily  that  D  is  constant  also.  Therefore 
the  likelihood  is  constant  for  all  y  >  0  .  As  an  illustration  of  the 
applicability  of  Theorem  5  we  consider 


£n\x  + 


U)'(x) 


_  1 
2  2 
(l+xZ) 


true  since  for  all  real  u 

3  5 

sinh(u)  =  u  +  +  ~y  + .  .  ■ 

Thus  we  can  conclude,  by  using  Theorem  5,  that 

u)(x)  =  4.  X(x)  =  P.n(x  +  J x2+l  ) 

satisfies  properties  (i),  (ii)  and  (iii). 

Notice  that  we  have  proved  that  the  likelihood  can  be  expressed  as 
the  sum  of  two  functions  which  are  of  the  following  form: 


Can  we  conclude  that  L  =  C  +  D  has  at  most  one  maximum?  We  cannot, 
because  we  must  assert  that  L'  =  C'  +  D'  has  at  most  one  zero  using  only 
the  fact  that  C'  >  0  and  D'  <  0  .  This  conclusion,  of  course,  is  impos¬ 
sible  based  only  upon  these  assumptions  (counter  examples  can  be  easily 
given) . 
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What  the  picture  suggests  is  that  C  and  D  must  be  concave  and 
convex  respectively.  If  this  were  true  it  would  be  more  likely  that  at 
most  one  maximum  could  occur.  An  effort  to  find  conditions  which  are 
sufficient  to  insure  those  properties  for  C  and  D  is  made  in  the 
next  section. 
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SECTION  VIII 


FURTHER  EFFORTS  TO  IDENTIFY  THOSE  DISTRIBUTIONS  FOR  WHICH  THE 
MAXIMUM  LIKELIHOOD  ESTIMATORS  DO  NOT  EXIST 

Before  we  proceed  we  collect  some  results  which  will  be  used  sub¬ 
sequently. 

Lemma  1:  If  <p  is  convex  increasing  and  H  is  convex  then  4>(H)  is 
convex  and  4(H)  increases  or  decreases  accordingly  as  H 
increases  or  decreases. 

Proof:  We  examine  the  first  t  /o  derivatives  of  the  composite  function, 

namely  .. 

4)'  (H)H’  ,  H"4>'(H)  +  (H')V'(H)  (34) 

and  the  result  follows  since  4'  ?  0  ,  <p"  >  0  and  H"  >  0  .  j  | 

We  have  the  immediate 

Corollary  1:  If  4  is  concave  increasing  and  H  concave  then  <J>(H)  is 

ccncave  and  decreases  or  increases  accordingly  as  H  decreases 
or  increases. 

Proof:  by  assumption  <}>'  >  0  ,  <p"  5  0  ,  H'  <  0  and  the  results  follow.  |[ 

Lemma  2:  If  <f>  is  convex  decreasing  and  H  concave  increasing  then  4(H) 
is  convex  decreasing. 

Proof:  Since  ty'  <■  0  ,  <p”  ?  0  ,  H'  >  0  ,  H"  n  0  ,  we  merely  examine 

the  expressions  in  (34). 
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We  now  examine  several  different  sets  of  assumptions  concerning  the 
behavior  of  (o  which  show  that  the  likelihood  function  when  y  is  fixed 
may  not  have  a  maximum  as  a  function  of  y  >  0  .  That  is  to  say  the  maxi¬ 
mum  likelihood  estimates  may  not  exist  for  all  three  parameters  jointly. 
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Let  us  consider  the  two  assumptions: 

(iv)  w' (x)  is  concave  for  x  >  0 

(v)  uj  *  (1/x)  is  log-concave  for  x  >  0  . 


One  checks  that  neither  of  these  assumptions  implies  the  other. 


From  assumption  (i)  we  know  that  d> (x ) 


ulSsI 


is  decreasing  on  (O,00) 


since  xcf) 1  ■  lj  -  —  <  0  .  By  assumption  (ii)  we  know  u)  is  concave  on 
(0,°°)  ,  and  thus  by  Corollary  2  it  follows  that  for  each  i  *  l,...n  the 
function 


h 


,  (y)  =  y .  |  4> ( I  y  I  / Y )  is  concave  increasing. 

1  lx 


Moreover  since  h^  >  0  we  know  is  also  concave  increasing.  Thus 


i  X  hi« 


is  concave  increasing. 


Now  since  -  —  £n  x  is  convex  decreasing  we  conclude  by  Lemma  2  that 


A(Y)  =  -  j  An 


;  X  hiw 


is  convex  decreasing. 


(35) 


Moreover  we  see  it  has  limits  of 


A(°°)  =  -  ~  In 


*2 


A(0+)  *  A(°°)  -  £n  u>'(°°)  . 


By  assumption  (iv)  we  also  know  that  to'  (x)  is  concave  decreasing.  Thus 
by  Lemma  3,  H(y)  =  aj'(|y^|/y)  is  concave  increasing.  But  £n  x  is  con¬ 
cave  increasing,  hence  by  Corollary  1  it  follows  that 


B(y)  =  ±  ^  An[w,(|y1|/y)] 


(36) 


is  concave  and  increases  to  zero. 
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Now  we  state 

Theorem  11:  Under  assumptions  (i) ,  (ii),  (iii)  and  (iv)  the  likelihood 
function  defined  in  (2),  with  fixed,  is  of  the  form 

L(Y)  -  A(y)  +  B(y)  for  y  >  0 

where  A  ,  as  defined  in  (35),  is  convex  decreasing  to  a 
finite  negative  limit  and  B  ,  as  defined  in  (36),  is  concave 
increasing  to  zero  and  moreover  L  has  no  maximum. 

We  defer  the  proof  until  later  and  state 
Theorem  12:  If  assumption  (iv)  is  replaced  by  (v)  in  Theorem  11  the  conclu¬ 
sion  remains  the  same. 

Proof:  We  have  only  to  check  that  (v)  implies  that  the  functions  defined 
for  y  >  0  ,  for  i  ■  l,...,n  by 

£n  w'dyJ/Y) 

are  concave  increasing  to  zero  and  hence  B(y)  is  also  concave  increasing 
to  zero.  |  ] 

We  now  give  the  additional  details  for  the  proof  of  Theorem  11.  Since 
L  *  A  +  B  with  A  as  defined  in  (35)  and  B  as  in  (36)  we  find 

L'  = 

Now  L'  >  0  iff 

n  n  n  ..u/i_.  i  /..\  /i__  i  /„.2 

^X 

J-l 


1  y  w 

n  La  i  i 


’j  nX 


"  V  I  /.(  I  '  I  /  \  I  I  /  I  / 


w'  (|y4 


Ehihi 

EhJ 


-iX 


^"( i y± I /y) ( |y± I /y2) 


“(|y. 
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"vtrrnim** 1  rr  •*«>  w- 


Multiplying  through  by  y 

and  letting 

*i  '  "l  ■  \ 

hi 

-•‘(lyJ/Y)  ! y±  1 

hi  * 

Ci  “  o),(|y1|/y'  y 

• 

we  can  write  (37)  as 

±  V  a 

n  L-* 

ihL  >  ^  1  ai 

•  i  V  «.  . 

n  . —  i 

(38) 

But  by  the  proof  of  Theorem  7  we  know  that  equation  (31)  follows  from 
assumptions  (i) ,  (ii)  and  (ill),  namely 


-  y  a  b.  >  -  \  a,  -N  b 
n  Z-  i  i  n  ZZ  in*_.  i 


We  can  establish  (38)  by  showing  that  b^ 


for  i  =  1 , .  .  .  n 


where  by  definition 


w"(|y. 1/Y)  |y. ! 

1  -  '  S'TTT'yT  ’  ~T 


This  last  inequality  is  true  if 


1  -  C(x)  >  -  ■-■77  T"  for  all  x  >  0  . 

u)  (x) 

but  since  by  (iii),  1  -  Q(x)  >  0  while  also  by  (ii)  co"(x)  <  0 
for  x  >  0.  we  have  the  last  inequality  established. 


Let  us  consider  again  the  case 


w(x)  =  £n(x  +  /x2+l )  . 


00  <  X  <  00 


Now 


n  (,)  -  li»  ^  -  y  li« 

U  co(x)  0)  (x) 


=  y  lim 


'1  +  x 


=  1 


x  -*■  00  /  1  +  (xy) 
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Thus  we  see  Che  function  D  is  negative  decreasing  with 


D(0+) 


On  the  other  hand 


6 


lim 

x  -*■  oo 


lira 
x  -*•  00 


X 


«*L 'Jsl 

0)(x) 


lim 
x  -*■  oo 


0  . 


As  a  consequence  the  negative  and  increasing  function  C  has  the  range  of 
values 

C(0+)  *  -  oo  f  C  (°°)  -  0  . 

Thus  L  =  C  +  D  and  C,  D  appear  schematically  as  drawn  in  Figure  3  . 
We  check  that  condition  (iv)  is  not  satisfied  for  all  x  >  0  since 
u>'"(x)  <  0  only  if  0  <  x  <  l//l~  .  To  see  this  examine  equation  (33)  . 

Now  we  check  to  see  if  condition  (v)  is  satisfied.  We  must  show  that 

Jin  oj'(1/x)  is  concave. 

Taking  the  second  derivative  and  substituting  from  (33)  we  find  that  it  is 
negative  iff 

-1  1-x2  2 

—  S  — — -r  iff  -  1  5  3x  for  x  >  0  . 

x  (1+x  ) 
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SECTION  IX 


CONCLUSION 

In  tills  report  we  have  examined  the  fatigue  process  and  broken  it  down 
into  three  distinct  intervals  each  with  separate  behavioral  properties. 

This  model  with  its  resuiting  complexity  has  been  adopted  as  an  alternative 
to  the  usual  stochastic  models  with  two-dimensional  parametric  distributions 
or  the  deterministic  models  involving  crack  growth  as  a  solution  of  a  differ¬ 
ential  equation  which  has  been  formulated  from  the  assumptions  of  fracture 
mechanics . 

A  three-dimensional  parametric  model  was  synthesized  which  incorporated, 
as  far  as  possible,  the  main  features  of  all  distinct  phases.  The  third 
parameter,  called  the  flexure  parameter,  is  an  alternative  to  introducing  a 
location  parameter  in  a  two-dimensional  model,  such  as  the  Wiebull,  before 
which  the  probability  of  failure  Is  zero. 

The  main  objection  to  the  introduction  of  a  location  parameter  as  "safe 
life"  is  that  if  such  a  time  does  not  exist,  in  fact,  then  serious  errors 
can  be  made  in  the  determination  of  the  initial  inspection  periods. 

The  three-dimensional  model  which  was  adopted  here  includes,  as  special 
cases,  many  of  the  simpler  two-dimensional  models  previously  studied.  Natur¬ 
ally,  such  a  general  class  must  be  carefully  substantiated  and  the  estimation 
techniques,  in  order  to  be  useful  in  the  specification  of  reliability  perfor¬ 
mance,  must  be  derived  for  the  types  of  fatigue  data  which  are  available. 

One  of  the  surprises  of  this  investigation  is  that  there  is  a  large 
class  of  such  distributions  which  conform  with  the  extant  data  for  which 
all  three  parameters  cannot  be  simultaneously  estimated  by  the  statistical 
method  of  maximum  likelihood. 
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The  study  of  this  last  section  was  completed  to  discover  those  conditions 
under  which  the  derivative  of  the  likelihood  function  for  the  unknown  flexure 
parameter  y  ,  when  the  other  two  parameters  are  unknown,  has  no  zero.  As  an 
important  practical  consequence  we  are  faced  with  the  situation  where  a  complete 
determination  of  the  parameters  determining  the  fatigue  process  cannot  be  accom¬ 
plished  by  maximum  likelihood  estimation  techniques.  Thus  certain  of  the  par¬ 
ameters  must  be  calculated  by  methods  which  rely  upon  knowledge  from  other 
disciplines.  Hence  those  disposable  constants  which  appear  in  the  model  must  be 
related  by  theory  to  constants  determined  from  the  material  or  the  geometry  of 
the  specimens. 
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LIST  OF  SYMBOLS 


Che  parameter  called  characteristic  fatigue  life 

the  regression  equation  for  characteristic  life  as  a  function 
of  maximum  stress :  the  Wohler  equation 

the  random  number  of  cycles  until  fatigue  failure 

the  total  random  crack  length  at  the  nth  cycle 

the  random  incremental  crack  growth  for  ith  cycle 

the  mean  and  variance  of 

the  standard  normal  distribution 


(with  or  without  affixes)  generic  symbols  for  distributions 


density  corresponding  to  upper  case  distribution 


shape  parameter 

(with  or  without  affixes)  the  critical  number  of  cycles  until 
failure 


the  disposable  function  expressing  the  reciprocal  relation¬ 
ship 

(witli  or  without  affixes)  the  random  time  for  a  phase  of  crack 
growth 


the  stochastic  crack  length  at  time  t 

the  ^-normal  variate  expressing  the  rate  of  crack  growth 


p(t  :  a,b)  the  propagation  factor  determined  as  the  solution  of  a  dif¬ 
ferential  equation  from  principles  of  fracture  mechanics 


A  the  ^-normal  variate  controlling  the  random  time  in  the 

third  phase  of  crack  growth 


the  standard  normal  variate 


<P-S 


-1 


M.Y.a 


the  location,  scale  and  flexure  parameter  for  the  log-life 
model 


w(x)  «  5(e)  the  disposable  function  for  the  log-life  variates 


the  likelihood  function 


P,R  product  and  ratio  functions  of  certain  derivatives  of  u) 

C,D  and  decompositions  of  the  likelihood  function 

A.B 

£  product  of  certain  derivation  of  to  indicative  of  the  degree 

of  s-shapedness 


h 


a  function  defined  using  a  reciprocal  of  the  argument  in 
u)(x)/x 
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